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Abstract 

We carry out driven diffusion Monte Carlo simulations of the two dimensional 
classical lattice Coulomb gas in an applied uniform electric field, as a model 
for vortex motion due to an applied d.c. current in a periodic superconducting 
network. A finite-size version of dynamic scaling is used to extract the dy- 
namic critical exponent z, and infer the non-linear response at the transition 
temperature. We consider the Coloumb gases / = 0, and / = 1/2, corre- 
sponding to a superconducting network with an applied transverse magnetic 
field of zero, and one half flux quantum per unit cell, respectively. 
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I. INTRODUCTION 



Phase transitions in two dimensional {2d) superconducting networks, such as peri- 
odic Josephson junction arrays and superconducting wire nets, has been a topic of much 
investigation.0 Theoretically, the phase transitions in such systems have been most exten- 
sively studied by equilibrium calculations and simulations.i^ Experimentally however, it 
has been most common to measure steady-state current- volt age (/ — V) characteristics, and 
look for a cross-over from linear to non-linear resistivity as a signal of the superconducting 
transition.i§ In this regard, the Kosterlitz-Thouless (KT) mo del§ of a vortex pair unbinding 
transition makes a clear prediction:0'lli as one heats up through the transition temperature 
Tkt, the I — V characteristics should make a discontinuous change from V ^ P exactly 
at Tkt to V ~ / above T^t- In experimental studies of superconducting 2d and 
networks! however, as well as in numerical simulationsjl^ agreement with this prediction has 
been claimed in some cases, not found in others, and is ambiguous in yet others. In particu- 
lar, it is not clear how this prediction may become modified when a transverse magnetic field 
is applied to the sample. In this case, the melting of the ground state vortex lattice induced 
by the magnetic field, may change the universality class of the superconducting transition, 
and lead to different steady state behavior. In view of these questions, it remains of interest 
to establish what steady state I — V behavior may be expected at criticality, for specific 
simple cases. 

Recently, a new dynamic scaling conjecture was proposed by Fisher et al^ to describe 
I — V characteristics in superconducting systems. Although this approach was developed for 
application to the "vortex-glass" transition in high temperature superconductors, it should 
apply equally well to any superconducting transition which is believed to be second order. 
In this work we carry out steady-state "driven diffusion" Monte Carlo simulation^! of the 
2d lattice Coulomb gasi'i in a uniform applied electric field, as a model for vortex motion 
due to a uniform applied d.c. current, in a periodic superconducting network. We consider 
the special cases of / = 0, corresponding to no transverse magnetic field, and / = 1/2, 
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corresponding to a transverse magnetic field of one fialf flux quantum per unit cell of the 
network. We apply a finite-size version of the new dynamic scaling conjecture to analyze 
our data, and extract the dynamic critical exponent z, and the power law of the I — V 
characteristic at the superconducting transition. For / = 0, we find z = 2, consistent with 
the Kosterlitz-Thouless prediction. For / = 1/2, we again find z = 2, consistent with the 
KT model, but in disagreement with expectations from recent equilibrium simulations of 
this model. 

The remainder of this paper is organized as follows. In Section II we give the theoretical 
framework for our simulations. We present our Coulomb gas model and the driven diffusion 
Monte Carlo algorithm. We review the KT vortex pair unbinding prediction, and we describe 
the finite-size version of dynamic scaling. In Section III we present our numerical results for 
/ = and for / = 1/2. In Section IV we give our conclusions. 



The standard modeloQ to describe behavior in a 2d superconducting network, is the 
uniformly frustrated 2d XY model, given by the Hamiltonian, 



Here 6i is the phase of the superconducting wavefunction at node i of the periodic network, 
the sum is over all nearest neighbor bonds (ij) of the network, and 



is the line integral of the magnetic vector potential A across the bond from i to j. For a 
uniform magnetic field applied transverse to the plane of the network, the sum of the A^j 
around (going counter-clockwise) any unit cell is. 



II. THEORETICAL FRAMEWORK 



A. The Driven Diffusive Lattice Coulomb Gas 




(1) 




(2) 
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J2A,, = 27Tf, f = BA/%, (3) 

cell 

where A is the area of a unit cell, and the constant / is the density of magnetic flux 
quanta ($o = '^e/hc) per unit cell. / is referred to as the "uniform frustration". U{(f)) 
is the interaction potential between the nodes of the network, which is periodic in with 
period 2-71, and has a minimum at = 0. For a Josephson junction array, one take^3 
f/(0) = — Jo cos(0). For a wire net, the Villain,lli or "periodic Gaussian" interaction may be 
more appropriate.^! 

It is generally believed that it is the excitation of vortices in the superconducting phase 
6i, that is responsible for the superconducting transition in such networks.l Since 2d vortices 
interact with a logarithmic potential,^ the Hamiltonain (|l|) is assumed to be in the same 
universality class (for the Villain interaction,lii the mapping via duality, is exact0) as the 
following Hamiltonian for Coulomb interacting point charges, 

'HcG = \ T.i^^ - f)V{n - - /)• (4) 

Here, i and j label the dual sites of the original superconducting network, = 0, ±1, ±2, ... 
is the integer vorticity of the superconducting phase at site and V{r) is the 2d lattice 
Green's function, which satisfies, 

AV(rl - r,) = -2Ti5i, (5) 

where is the discrete Laplacian. In this work we restrict our interest to a square network, 
with periodic boundary conditions. In this case, V{r) is explicitly given by 

v{f) = ^ y Z^-^' ^ (6) 

^ ' N ^ 2 - cos - cos ky ^ ^ 

where k are the allowed wave vectors with k^ = {2'Kn^/L), with = 0, 1, L — 1. L is the 
length of the network, and N = L"^ . For large r, 

V{f)^\nr. (7) 

Since V{f = 0) is divergent, the partition sum over {nj} is restricted to neutral configurations 
where 



^E-^ = /- (8) 

i 

See Ref. ^ for further details. Thus the average density of vortices is equal to the density of 
flux quantum of the applied magnetic field. 

The Hamiltonian @) therefore represents a density / of integer point charges, on a 
uniform compensating background charge — /, interacting with the 2d Coulomb potential. 
For / = 0, the ground state is the vacuum, and low lying excitations are bound neutral pairs 
of rii = ±1. For / = 1/2, the ground state is an ordered checkerboard lattice of chargeJl 
with rii = +1 on every other site, as shown in Fig. 1. Low lying excitations may be viewed 
either as a displacement of one of the charges in the ground state lattice, or equivalently 
the creation of a neutral Arij = ±1 pair of charges superimposed on the ground state lattice 
configuration. It is this Coulomb gas analog of the superconducting network which we will 
use to carry out our simulations. 

To model flux flow resistance in the superconducting network, we apply a uniform electric 
field to the Coulomb gas charges nf, this models the uniform Magnus force0 that an applied 
d.c. current exerts on vortices in the superconductor. The net charge current density in 
the Coulomb gas then corresponds to the net vortex current density in the superconductor, 
which is proportional to the net flux flow electric field (or voltage drop per unit length) in 
the superconductor. If we denote as E the applied electric field, and J the resulting charge 
current density, in the Coulomb gas analog, and JT" the applied d.c. current density, and £ 
the resulting flux flow electric fleld, in the superconducting network, then the correspondence 
between the models is given by 

E^J, J^£. (9) 

The simulation of the Coulomb gas in the presence of the uniform E is carried out by use 
of the "driven diffusion" technique.0 We imagine adding to the Hamiltonian (^) the term 

Sn = -Y.nifi-E (10) 
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which represents a dipole interaction between the charges and the electric field E. Al- 
though HcG + is unbounded, and therefore not a valid Hamiltonian in a global sense, 
when it is used as a local energy function for computing energy differences, in connection 
with the standard Metropolis Monte Carlo algorithm for accepting or rejecting proposed 
excitations, it yields an enhanced probability (consistent with detailed balance) for accept- 
ing excitations with a net movement of charge in the direction of E, thus setting up a 
non-equilibrium steady-state with a finite charge current density J fiowing parallel to E. 

In our simulations, we have chosen E = Ey, along one of the axes of the periodic lattice 
of sites. At each step of the simulation we pick at random a pair of nearest neighbor sites 
{io, ii). For the / = case (where the ground state is rii = 0), we add a positive unit charge 
to site iq, ie. Ariig = +1, and a negative unit charge to site ii, ie. Arij^ = —1. For the 
/ = 1/2 case (where the ground state is as in Fig. 1), we simply interchange the charges, riig 
and njj, at the two sites. For / = 1/2, this restricts configurations to those where half of 
the Hi are +1 and the other half are 0; charges Ui = —1 or +2 are not allowed. Tests showed 
that these other values of Ui correspond to higher energy excitations, which are negligible at 
the temperatures of interest, ie. the melting temperature of the ground state charge lattice. 
In both the / = and the / = 1/2 cases, the change in energy due to the addition of the 
excitation is computed using Hcg + ^'H, and the excitation is accepted or rejected using 
the Metropolis algorithm. In both cases, acceptance of the excitation gives a contribution 
to the average current density, 

AJ = An J''~'''' + AnJ-^^^ = An,,if,, - rlj (11) 

where Aui is the change in rii at site i created by adding the excitation, and our algorithm 
always satisfies Arii^ = —Arii^. 

While the above driven diffusion Monte Carlo algorithm encodes a specific dynamics, 
which in detail may well be different from the true microscopic dynamics of vortices in 
superconductors, our hope is that qualitative behaviors which are largely determined by 
energetics, particularly the non-linear form at criticality, will be preserved. We have chosen0 
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to simulate the driven diffusion Coulomb gas, rather than the more realistic RS J dynamics for 
an XY model of a Josephson array, because in the Coulomb gas algorithm one directly 
moves the important degrees of freedom, the positions of the vortices. This results in a 
computationally faster algorithm for several reasons: (z) spin wave degrees of freedom are 
eliminated; (ii) the effective energy barrier for an isolated vortex to hop to a neighboring cell 
in the XY mo del0 is removed, since in the Coulomb gas this hop takes place in one discrete 
step; (iii) the RSJ dynamics requires a computation of order (or L\nL for improved 
algorithms^) at each step of simulation, where the Coulomb gas requires a computation of 
order only when the trial excitation is accepted;i for the low acceptance ratios we find, 
this effect is significant. 

N = L? steps of the above updating process will be referred to a one MC pass. In our 
runs, presented in Section III, an initial 10,000 passes were typically discarded to equilibrate 
the system. Following this equilibration, five independent runs (using different random 
number sequences) of 200,000 passes each, were used to compute averages. Our error bars 
are estimated from the standard deviation of the averages from these five runs. 



B. Kosterlitz-Thouless Pair Unbinding Model 

We now review the Kosterlitz-Thouless model of pair unbinding, as applied to the 
determination of non-linear steady-state behavioi0 below the transition temperature Tkt- 
If we consider the addition of a neutral pair of charges An^ = ±1 separated by a distance r, 
we may estimate the free energy of this pair in the presence of all other charges as, 

Fpair{^ = -^{\n\r\-E-f). (12) 

Here e is the effective long wavelength dielectric function of the Coulomb gas, which serves 
to screen the logarithmic interaction between the members of the pair, as well as to screen 
the interaction of the pair with the external field E. A pair oriented along the direction of E 
therefore sees a potential maximum at = IZ-E. If the pair is able to overcome this potential 
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barrier through thermal fluctuations, the pair can then lower its free energy by increasing 
the separation r without bound. The pair will thus unbind, and give a contribution to a net 
current of charges flowing along the direction of E. The rate for such critical pair unbindings 
to occur is given by the Boltzmann factor, 

Wpair ~ e-^-"-(^^)/^ ~ E^/^^. (13) 

Such critical pairs will expand until they recombine with other such free charges into new 
bound pairs. This unbinding and recombination of pairs leadJi3 to an effective density of 
free charges Uf, 

Uf ~ {WpairY^^ (14) 

Using Eqs. ([13|Jl^) , the net current that flows due to pair unbinding is then. 

The insulator to metal transition in the Coulomb gas, where e diverges, marks the cross over 
from non-linear to linear J — E characteristics. Using the correspondence of Eq.(^, together 
with / = LJ and V = LS for the total current and total voltage drop in a superconducting 
network, we see that this Coulomb gas insulator to metal transition corresponds to the 
superconducting to normal transition in the superconducting network. 

For / = 0, where the ground state is the vacuum, pair unbinding excitations such as 
above, are believed to be the only source of net charge current. The Kosterlitz-Thouless 
model is expected to describe the insulator to metal transition in this system, and gives 
the prediction0Ei that e^^(T) has a universal discontinuous jump to zero exactly at the 
transition temperature T^t, with, 

l/e{TKT)TKT = 4. (16) 

Thus, exactly at T^t, Eq.(pr5D gives the non-linear behavior, J E^. The corresponding 
result in the superconducting network is V ^ . 
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For / = 1/2, where the ground state is the doubly degenerate ordered charge lattice 
shown in Fig. 1, the above pair unbinding continues to provide a mechanism for non-linear 
response below the insulator to metal transition temperature, which we will continue to refer 
to as Tkt- However there are now also other possible excitations, involving the formation 
of domains of oppositely oriented ground state, which might possibly contribute0 to a non- 
linear response in J. Thus no clear prediction exists for the form of the non-linear response 
at the transition. 

Similarly, the nature of the equilibrium transition in the / = 1/2 model remains 
controversiahlii'iS If Tkt is the insulator to metal transition, the Kosterlitz-Thouless argu- 
ments continue to provide a lower bound on a discontinuous jump in e~^, ie. 1 / e{TKT)TKT > 
4. It is unclear however, whether this is satisfied only as an inequality, or whether the uni- 
versal KT behavior as in Eq.([16|) continues to hold. Additionally, there is a well defined 
temperature Tm in the model, corresponding to the melting of the ordered ground state 
charge lattice. General arguments^ suggest the bound Tm > Tkt^ however it remains un- 
clear whether or not these two temperatures are in fact equal. It is further unclear whether 
the melting transition at Tm is Ising-like (as the double discrete symmetry of the ground 
state suggests), or whether the long range interactions between the charges cause the melt- 
ing to fall in a new universality class. Our present work was in part motivated by the hope 
that dynamic calculations might shed some light on these remaining equilibrium questions. 

C. Finite-Size Dynamic Scaling 

Recently, Fisher et a/.0 have proposed the following dynamic scaling relation, for an 
infinite superconducting system with a second order phase transition at Tc- The relation 
between the dissipative electric field and the applied d.c. current density is given by, 

£ = Ji''-^-'^^{JE,''-^/T) (17) 

where ^ is the spatial correlation length, d the dimensionality of the system, z the dynamic 
scaling exponent, and <l>± the scaling function above and below the transition. The most 
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natural generalization of this form, to a system of finite length L, is, 

S = Jh'^-^-'^Jh'^-^/T, tb^/", b/L) (18) 

where b is an arbitrary length rescaling factor, t = {T — Tc)/Tc, and u is the correlation 
length scaling exponent, ^ ~ t^". The form Eq.(|T7p can be obtained from Eq.(|IBp by 
choosing b = C,, and taking L — * oo. For finite L, $ is a continuous function as t passes 
through zero. Only in the L ^ oo limit does ^{J',t,0) become discontinuous as t passes 
through zero; this determines the different scaling functions and $_ of Eq. ([l7|) . 

From Eq.(|l8]) one can now determine the scaling behavior at criticality, t = 0. Choosing 
b = L — > oo, and setting t = 0, one finds, 

as found by Fisher et al^ Thus at Tc, one always expects a non-linear power law response. 
For d = 2, a. dynamic exponent of z = 2 recovers the S ~ J^^ result of the KT pair unbinding 
picture. 

Using the correspondences of Eq.(^, we now recast the scaling equation ( [T^ ) into a form 
for use with our Coulomb gas model. Choosing the rescaling length b = L, we get,ii 

J = EL'^-^-'^^EL'^-^/T, tL^'^ 1). (20) 

Finally, for d = 2, exactly at criticality, t = 0, we have the scaling, 

J = EL-'^{EL/n,0,l). (21) 

To fit this scaling equation to our Monte Carlo data, and extract the dynamic exponent z, 
we use the method used by Nightingale and Blote@ for similar equilibrium problems. We 
consider behavior exactly at Tc, as a function of varying E, in the limit of large L but small 
EL. Expanding the scaling function $ gives, 

J{E, L) = EL-'[(^o + ^lEL + ^2{ELy + 0{ELf]. (22) 
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Truncating this expansion at any finite order, we perform a least non-linear fittZl of our 
Monte Carlo data to Eq.(p2D to determine the unknown parameters 2;, and the 

We check for stability of our fit by increasing the order of the expansion, and by dropping 
data from successively smaller values of L, and checking if the fitted parameters change 
within our estimated statistical error. Statistical errors in the fitted parameters are estimated 
by generating many synthetic data sets, by adding random noise to each of the original MC 
data point. The noise for each data point is taken from a Gaussian distribution whose width 
is set equal to the estimated statistical error of the original MC data point. Using these 
fictitious data sets, we repeat the fitting procedure to obtain new fitted parameters. The 
estimated error of a parameter is then taken as the standard deviation of the results from 
all the fictitious data sets. 

III. NUMERICAL RESULTS 
A. / = 

For our simulations of the / = Coulomb gas, corresponding to the ordinary XY model, 
we use as the equilibrium KT transition temperature Tkt = 0.218, as determined by one 
of usS from a finite size scaling analysis applied to equilibrium simulations of e~^(T, L). 
This value is in good agreement with earlier estimates, based on Coulomb gas simulations 
by Saito and Miiller-Krumbhaar,i Tkt = 0.215, and by Grest,i Tkt = 0.220. In Fig. 2, 
we plot the resulting charge current density J{E, L) versus for several values of L, at 
the fixed temperature Tkt = 0.218. We see that the smaller E, the larger is the finite size 
effect as L varies. From Eq.(^I]) we see that smaller values of E probe larger length scales; 
correspondingly, our statistical error increases as E decreases. 

To find the dynamic exponent z, we now fit the data of Fig. 2 to the expanded scaling 
function of Eq. (]22|) . Since this expansion converges fastest for small values of the argument, 
we restrict the data used in our fitting to those points where EL < 1. This corresponds to 
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lattice sizes L = 6 — 14, with E = 0.02 — 0.08. In Table I, we show the results of this fit, 
for several orders of expansion, for various ranges of L. Using the fourth order expansion 
for lattice sizes L = 8 — 14, we find z = 2.073 ± 0.098. Using this z in Eg. (pj]) , we get a 
non- linear response J ~ ^3.073^ consistent with the prediction from the KT pair unbinding 
model, assuming the universal jump in e~^(TxT) (see Eqs.(^,|I^)). 

In Fig. 3 we plot our data as JL^/E versus EL. We see that the data collapses onto a 
universal curve representing the scaling function 0, 1). The value z = 2.073, obtained 
from the fitting, is used in making the vertical axis, and the solid line is drawn using the fitted 
values of the $j. Although the agreement is reasonable. Table I does suggest some potential 
problems. The parameters 2; ~ 2 and $0) while remaining stable within the estimated 
errors, both show a systematic increase as the smallest Lj used in the fit is increased. $2, 
although consistent within the different fits, shows a very large estimated statistical error 
(other $j, although also strongly fluctuating, seem too small to be significant in the fit). 
Ideally, one would like to carry out these fits using increasingly smaller values of E than we 
have used here. However, when E becomes small, equilibration times become large, and we 
were unable to get accurate enough data to improve our fit. 

Our results in Table I and Fig. 3 represent checks of scaling in the small EL limit. We 
have also tried to check scaling in the large EL limit. Provided that EL is sufficiently 
large that finite size effects are small, we should be able to use Eq . ([ITD to collapse our data 
onto two universal curves, given by and $_ above and below T^t, by plotting J^^ /E 
versus E^/T. To do so, we need an expression for the correlation length ^(T). For the 
Kosterlitz-Thouless transition, asymptotically close to T^t^ this form is,01li 

e±(T) ~ e^±/l^-^^-l^ (23) 

where the subscripts + and — refer to behavior above and below T^t respectively, and for 
the KT transition, z/ = 1/2. Minnhagen and Olssonil have argued that this true asymptotic 
form holds only in a narrow critical region of about 5% of Tkt- Nevertheless, they also 
indicate that Eq. (^) is a useful phenomenological form for fitting over a wider temperature 
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range for T > Tkt, provided C+ is taken as a phenomenological parameter not necessarily 
equal to the true asymptotic value. We adopt this approach and use Eq.(p3|) with C± and 
u as phenomenological parameters. 

To carry out this large EL check of scaling, we observe from our data at Tkt in Fig. 
2, that finite size effects are negligible provided we restrict the data to L > 24, E > 0.06. 
Since this is true at Tkt, it should also certainly be true for other values of T. We therefore 
carry out simulations on an L = 24 lattice, for values E > 0.14. Our results for J{E,T) 
versus E, for various T above and below Tkt, are shown in Fig. 4. on a log-log scale. Solid 
lines with slopes of 1 (for Ohmic behavior above Tkt), and 3 (for critical behavior at Tkt) 
are shown for reference. In Fig. 5, we try to collapse this data onto two universal curves 
as discussed above, by finding the best choices for the parameters Tkt, z, C±, and z/. The 
results shown are for the values Tkt = 0.218, z = 2 (consistent with our small EL analysis), 
u = 1/2 (consistent with the KT form), and C+ = C_ = 0.35. The collapse is reasonable, 
except for the smallest several values of T below Tkt- This is most likely due to a failure of 
the assumed form for ^(T), Eq. (P5D , to be valid over such a large temperature range. 



B. / = 1/2 

In this section we carry out a similar analysis as in the previous section, except applied 
to the / = 1/2 Coulomb gas, which corresponds to the "fully frustrated" XY model. As 
discussed at the end of Section ^IB[ there are in principle two transitions in this model: a 



insulator to metal transition at Tkt, and a charge lattice melting transition at Tm- It is 
Tkt which corresponds to the transition from non-linear to linear J — E characteristics (see 
Eq. (pISD and following discussion). The most recent equilibrium simulations of the / = 1/2 
Coulomb gas model by one of u^l find that Tkt — 0.126 is very close to, but slightly 
below Tm — 0.1315; the discontinuous jump in is l/e{TKT)TKT — 5.35, larger than the 
universal KT value of 4 (see Eq. (p^) ) . This compares with earlier estimates by Gresti of 
Tkt = 0.129 ± 0.002, with a jump l/t{TKT)TKT ^ 4.88 ± 0.31. Similar simulations on 
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the fully frustrated XY model by Ramirez-Santiago and JoseQ find Tkt ~ Tm and a jump 
l/e{TKT)TKT ^ 5.21. 

Fixing the temperature at Tkt = 0.126, we show in Fig. 6 our results for J{E, L) versus 
E for various L. To extract the critical exponent z we fit this data to the expanded scaling 
function of Eq.(^). We restrict!^ the data used in our fitting to those points where EL < 0.5, 
which corresponds to lattice sizes L = 6 — 14, with E = 0.01 — 0.04. In Table II, we show 
the results of this fit. Using the fourth order expansion for lattice sizes L = 8 — 14, we find 
z = 2.060 ± 0.124. As was found for / = 0, the result z ~ 2 is consistent with a power law 
response of J ^ E^. 

In Fig. 7 we plot the data of Fig. 6 as JL^/E versus EL, and find fair agreement 
with the expected collapse onto a universal curve. Our fitted value of z = 2.060 is used 
in making the vertical axis, and the solid line is drawn using the fitted values of the 
Although the agreement is reasonable. Table II again suggests some potential problems. As 
we found in Table I for the / = case, now for / = 1/2, the parameters z and $o show 
a systematic increase as the smallest Lj used in the fit is increased. Now however, this 
increase is more pronounced, and the fitted values remain consistent with varying Lj only 
within the outer limits of the estimated statistical errors. Furthermore, the higher $j all 
seem to be significant, and all have very large statistical error. These observations make it 
unclear whether or not our data truly represents the asymptotic scaling region of large L, 
small EL. 

We have not attempted to check scaling for this / = 1/2 model in the large EL limit. 
The strong finite size effects seen in Fig. 6, even comparing L = 24 and 32, means that we 
would have to go either to larger lattice sizes (which are beyond our present computational 
ability), or to temperatures sufficiently far from Tkt, in order to reach the large EL limit, 
for the values of E we have studied. Uncertainties in the correct form one should take for 
C,{T), due in particular to the close proximity of the vortex lattice melting transition at Tm 
to the insulator to metal transition at Tkt, would undoubtedly make such an analysis more 
complicated than was the case for / = 0. 
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IV. CONCLUSIONS 



To conclude, we have carried out steady-state driven diffusion Monte Carlo simulations of 
the 2d lattice Coulomb gas, in order to compute the dynamic exponent z, and hence obtain 
the non-linear response J^E°-,a = z + l,aX criticality. This corresponds to the non-linear 
current- volt age characteristic \^ ~ J" in a superconducting network at the superconducting 
to normal transition. We have analyzed our data according to a finite-size scaling method 
based on a new dynamic scaling conjecture by Fisher et al^ 

For the / = model, corresponding to a superconducting network in zero applied mag- 
netic field, our results agree with the familiar Kosterlitz-Thouless pair unbinding model. Our 
finite-size scaling analysis, varying L and E at fixed T = T^Ti gives a value of z ^ 2, consis- 
tent with a power law response at Tkt of a = 3. Our check of scaling in the infinite L limit, 
where we vary T and E for fixed L large, shows fair agreement with the Kosterlitz-Thouless 
model, but success is limited by our limited knowledge of the form of the correlation length 
^(T) outside the narrow critical region. 

For the / = 1/2 model, corresponding to a superconducting network in an applied mag- 
netic field of one half flux quantum per unit cell, we again find z ^ 2. This is consistent 
with the Kosterlitz-Thouless pair unbinding result of Eg. (p!5D only if the discontinuous jump 
in e'^iTxT) obeys the universal KT prediction of Eq.(^), ie. l/e(TKT)TKT = 4. Equi- 
librium simulations however indicate that for / = 1/2, this jump is non-universal, with 
l/e{TKT)TKT — 5.35. Using this value of e in Eq.(|l^) yields the non-linear response at Tkt 
due to pair unbinding as, J ~ E°-, with a = 3.68. If pair unbinding were the dominant 
contribution to J, this would imply a dynamic exponent of z = a — 1 = 2.68. 



It remains unclear what is the source of this inconsistency. It could be that the analysid-l 
of in the equilibrium simulations is in error; this equilibrium analysis involves identifying 
leading logarithmic corrections at Tkt which may be hard to determine accurately. Or it 
could be that our finite-size analysis of z in this present work is flawed; this possibility is 
indicated by the less than satisfactory behavior of the fitted parameters (see Table II, and 
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discussion at the end of Section [III B| ) as we vary the order of the fitting expansion, or the 



system sizes used in the fit. Or it could be that the result 2; = 2 is a more general property of 
such superconducting systems, which is independent of the KT pair unbinding model; in this 
case one would expect that some excitation other than pairs gives the dominant contribution 
to J at Tkt- The natural guess for these other excitations is the domain excitations of the 
ground state charge lattice. However this would seem unlikely, if the charge lattice melting 
transition Tm is distinctly higher than Tkt-, as equilibrium simulations suggest. Therefore, 
behavior in the / = 1/2 model remains an enigma, both from the equilibrium, and now from 
the steady-state dynamic point of view. 
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FIGURES 

FIG. 1. Ground state charge lattice for the / = 1/2 Coulomb gas. A + denotes the presence 
of a charge iii = +1. 

FIG. 2. Charge current density J{E, L) versus applied electric field E for various square lattice 
sizes L, at fixed temperature Tkt = 0.218, for the / = model. 10^ total MC passes were used to 
compute averages. 

FIG. 3. The finite size scaling behavior of the charge current density. JL^/E versus EL is 
plotted for various system sizes L at fixed temperature Tkt = 0.218, for the / = model. Symbols 
with error bars represent the MC data. The solid line results from the fitting to Eq.(22) using a 
fourth order expansion in EL, and data from L = 8 — 14 and E = 0.02 — 0.08. Data from E = 0.1 
and 0.12 for each lattice size are included in the plot. The fitted value of z = 2.073 was used in 
making the vertical axis. 10^ total MC passes were used to compute averages. 

FIG. 4. Charge current density J{E, T) versus E for various temperatures T, for fixed system 
size L = 24, for the / = model. The solid lines from left to right have slopes 1 and 3 respectively, 
illustrating Ohmic behavior above Tkt, and KT critical behavior at Tkt- 10^ total MC passes 
were used to compute averages. 

FIG. 5. Test of the infinite size scaling relation Eq.(17) for the / = model. The data of Fig. 
4 is replotted as J^^/E versus E^/T, using the form of Eq.(23) to determine ^{T). Parameters z, 
Tkt, I', and C± are varied, to get the best collapse of the data onto two distinct scaling functions 
$-1-, above and below Tkt- 

FIG. 6. Charge current density J{E, L) versus applied electric field E for various square lattice 
sizes L, at fixed temperature Tkt = 0.126, for the f = 1/2 model. 10^ total MC passes were used 
to compute averages. 
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FIG. 7. The finite size scaling behavior of the charge current density. JL^ /E versus EL is 
plotted for various system sizes L at fixed temperature Tkt = 0.126, for the / = 1/2 model. 
Symbols with error bars represent the MC data. The solid line results from the fitting to Eq.(22) 
using a fourth order expansion in EL, and data from L = 8 — 14 and E = 0.01 — 0.04. Data from 
E = 0.05 and 0.06 for each lattice size are included in the plot. The fitted value of z = 2.060 was 
used in making the vertical axis. 10^ total MC passes were used to compute averages. 
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TABLES 

TABLE I. Results of the fitting of J to an expansion of the scaling function in powers of EL, 
as in Eq.(22). The data of Fig. 2 for the / = model is used. The first column shows the range of 
system sizes Li — Lf which are included in the particular fit. The following columns give the fitted 
parameters. The last column is the error of the fit. For each sequence of L^ — Lf, the first row 
gives the value of the fitted parameter, while the second row gives the estimated error. 
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TABLE II. Results of the fitting of J to an expansion of the scahng function in powers of EL, 
as in Eq.(22). The data of Fig. 6 for the f = 1/2 model is used. The first column shows the range 
of system sizes Li — Lf which are included in the particular fit. The following columns give the 
fitted parameters. The last column is the error of the fit. For each sequence oi L^ — Lf, the 

first row givew the value of the titled ]")arMUieler. while the scx'oud row gives the estiiualcd enx^r. 



Li - Lf 


z 




$1 


$2 


$3 


$4 


■N/2 

Xfit 








1st order 


expansion 






u — 1^ 


1 SQ7fi 


n 7788 


1.3337 










1 4 9fifi1 




u.uooy 


n 1 1 '^4 


0.6810 












8 — 14 
o — 1^ 


l.iJOOU 


n 88Q1 


1.5631 










1 9 ^81^9 






u. iooy 


0.8506 












in 1/1 


1 QSQ/I 


n Q^97 


1.7894 










1 n i^Qsn 




U.Uoyi 


U.UOM 


0.6560 


















2nd order 


expansion 






fi 1/1 


1 Q1 

i.yioo 


n QnQ9 
u.yuyz 


0.6484 




1.1746 






1 1 9ni^i 




U.Uc7Z J. 


n 1 7nfi 

U. -L I UU 


0.6392 




0.6416 








8 1/t 




1 4c;9'^ 


-0.2034 




3.4457 






4 p;49Q 




n 1 9^"^ 

U. iZOO 


n i^^fi/i 

U.0OD41: 


0.5575 




2.7420 








10 — 14 


9 Til4 


1 8496 


-0.5706 




4.8902 






3 0757 




0.0680 


0.2814 


0.9038 




1.7059 














3rd order 


expansion 






6-14 


1.9128 


0.8566 


1.2388 




-0.9075 


2.0827 




11.0736 




0.0938 


0.2125 


1.0552 




3.6274 


3.5002 






8-14 


2.0602 


1.3639 


0.7936 




0.0519 


3.3454 




4.3721 




0.1246 


0.3867 


1.8768 




4.4841 


6.9044 






10-14 


2.1581 


1.6954 


1.5878 




-2.4100 


7.3724 




2.5403 



24 





0.0699 


0.2873 


1.0701 


1.8182 


5.3878 






4th order expansion 


6- 


14 


1.9125 


0.8636 


1.0860 


0.0985 


-0.4950 


2.1906 


11.0117 






0.0896 


0.1444 


0.9954 


1.7584 


2.9752 


3.0877 




8 - 


14 


2.0597 


1.3482 


0.8797 


0.3090 


1.4246 


2.2639 


4.3701 






0.1238 


0.4363 


0.9554 


1.1924 


2.0356 


5.9653 




10- 


14 


2.1554 


1.6883 


1.3351 


-0.2117 


1.0093 


5.6638 


2.5421 






0.0821 


0.3368 


1.3246 


1.9165 


1.6536 


6.1369 





25 



